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SUMMARY 

This article summarizes some aspects of research in progress for develop- 
ing finite element methods for contact problems. We propose a new"finite ele- 
ment approach" for contact problems in two-dimensional elastodynamics . Sticking, 
sliding and frictional contact can be taken into account. The method consists 
of a modification of the shape functions, in the contact region, in order to 
involve the nodes of the contacting body. The formulation is symmetric (both 
bodies are contactors and targets), in order to avoid interpenetration. Compati- 
bility over the interfaces is satisfied. The method is applied to the impact of 
a block on a rigid target. The formulation can be applied to fluid-structure 
interaction, and to problems involving material nonlinearity. The extension to 
three dimensions presents additional difficulties, but it is possible. 


INTRODUCTION 


The approach presented in this article was developed while trying to simu- 
late the movement of a gas bubble in a liquid. The original idea was to intro- 
duce the compatibility of the velocities over the gas-liquid interface via a 
constraint equation and to handle it by the Lagrange multiplier method. In a 
second step, the Lagrange multiplier method was replaced by a penalty method, 
which is easier to implement. In both cases, the constraint equation is a 
geometric relationship between gas and liquid velocities. No local remeshing 
was performed; the bubble and liquid meshes were simply superposed. This resul- 
ted in poor pressure fields along the interface. Looking for an improvement of 
this situation, remeshing appeared as the best but also the most cumbersome 
solution. Alternatively, a modification of the shape functions appeared to have 
the advantages of remeshing, without its inconveniences. This latter approach 
is described herein as it is applied to contact problems in two-dimensional 
elastodynamics. Frictional contact results in an exchange of momentum between 
the two contacting bodies, and can be realised by direct introduction of a 
contribution of the contactor's velocity into the target's equation of motion. 
This is conveniently done by means of a modification of the shape functions, as 
described in the next paragraphs . 
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The proposed approach has the advantage , as compared to the Lagrange mul- 
tiplier method, of maintaining a constant size of the linear system to he 
solved. Compared to a penalty method, it has the advantage that we get auto- 
matic compatibility of the field variables over the interface. When the 
formulation is symmetric (i.e., both bodies are targets and contactors), 
interpenetration is totally avoided. 


MODIFIED SHAPE FUNCTIONS 
FOR QUADRILATERAL FINITE ELEMENTS 


Figure 1 shows a two-dimensional contact problem. Node C contacts element 
(1-2-3-4) of the target and from then on contributes to its shape functions . 

We start from the initial 4-nodes interpolation function 

V = i N|*, (!) 

a = 1 

with N a = 0.25 (1+ sign(g a ) g ) ( 1+ sign (7j a ) 7) ) (2) 

where S', 7] are local coordinates, and v are the velocities. 

In order to take the contribution of point C (node 5) into account, we modify 
the interpolation function as follows : 


J V 
a v a 


V = X N 3 
a=l 

Notice that from a global point of view there is no new node appearing. 


(3) 


Obviously, a "hat shape function" at node 5 is the most adequate for our 
purpose. This yields automatic compatibility of the velocities at the interface, 
if a symmetric formulation is used. Further, we want to account for tangential 
sliding with friction at the contact point. Therefore, we introduce a factor p. 
which allows the shape function at node 5 to vary in amplitude between 0 and 1, 
which will lead to a partial exchange of momentum. The resulting shape functions 
are (assuming C is on side 7) = + l) : 


N 


* 

5 




0.5^(l + 7,)[U(g-g 5 )/(l + § 5 )] 
0.5[i (l + 7|)[l-(g-g 5 )/(l-E 5 )] 


if g<E 5 
if 1>U 


(4) 


and 

with 


n; = N a - a, N 


(5) 


a. 


= 0.5 (l+sign ^ g 5 ) 


The local coordinates of C are (§ ,l). We can also assume, without restric- 
tion, that the contact point is associated with local coordinates (0, +l), N* 
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N* = 0.5 p, (l + f])(l - sign (E ) • l) 


( 6 ) 


then becomes 


t 


This shape function is shown on figure 2. Observe that N * (g,., 71^ ) hoes not 
vanish at node 5 when p,^ 1, but N = 1 is preserved. 

£L 

In the sequel, we separate normal (n) and tangential (t) directions and 
therefrom the following typical possible situations. 


a. p. = 1, (1 = 0, corresponds to frictionless sliding in the tangential direc- 

tion and sticking in the normal direction. This yields : 

* 

N = N (a = l-*-5 ). the standard 5~nodes interpolation 
n a 

function, 

% 

N = N (a = 1— U), the standard U-nodes interpolation 
9 » 

function. 


b. u = 1, p. = 1 corresponds to sticking, and yields : 
n L- 

N* = N* = N (a = 1— 5). 
n t a 


c. |l n — 1, p t e ]0.1 [ this accounts for frictional sliding. 

Since p depends on orientation, we introduce a second order tensor, which we 
need in order to define strain rates and stresses in global coordinates. 

In a local orthogonal frame tangential to the target surface we write 


*1 


M-n 

0 


v n C 



0 



V , C 


( 7 ) 


where the superscripts T and C stand for target and contactor, respectively, and 
the subscripts n and t for normal and tangential. 

, Equation (7) defines the contribution of node 5 (contactor) to the target 
velocity while the true local velocity, at C, is given by 


sign (g) 


-1 if g < 0 
+1 if g >0 
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)v n |=|K] M (8) 

where"*" v ° — [ v n v / J (9) 

I v n = I v T v T v T v T v T v C v C | T (10) 

| a [ In It v 2n 2t 4t 5n 5t J v ' 

and [ N " J is defined as follows"*"*" 

[N 5 n ]=N 5 [n n ] (11) 

K]= N a [i ] - a a [Ng] (12) 

M = [o" l\ <«) 


The matrix is diagonal when defined in a local referential, tangential 
to the contact surface. It does not induce a coupling of the normal and tangen- 
tial components, but this would not be true in a global referential. We can, 
therefore, establish the stiffness in this local referential and rotate the 
element matrix before we assemble the elements. Alternatively we can make the 
derivation in a global coordinate system. 

In practice, the whole effort essentially reduces to minor changes in the 
shape function routines. 

TRANSIENT SOLUTION PROCEDURE 
Search Algorithm 

We need to determine at each time step the location of each node of each 
body in the contact zone with respect to the mesh of the other one. For that 
purpose a connectivity matrix is established in the input phase; this matrix 
lists all elements connected to each element. Assuming the time step to be 
small, we memorize the previous position of each node (by an element number), 
and search for its new position in adjacent elements . A 2-dimensional search 
path is shown on figure 3. Once the new position is known, we modify the 
shape functions of the target element as described in the previous paragraph 
and compute the updated stiffnesses. 


t T as superscript of a matrix stands for transpose 

tt a depends on the local coordinate of the fifth node 
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Predictor-Corrector Algorithm 

We adopt here an explicit predictor-corrector algorithm, defined by the 
following equations, at time step (n + l) (see ref. 1, 2 for details). 


M J nt1 + N(d n *i »? n .l) = F n .1 

(15) 

d„.i — + At v n + —■ (l - 2 (3) a n 

( 16 ) 

?n.i = In + At(l-\)a n 

(17) 

£n.1 = £n*1 + A<2 P 3n»1 

( 18 ) 

Jin*1 = 5n.1 + At Ya n .i 

(19) 

d 0 =5 

(20) 

*0 =1 

(21) 

*0 = M- 1 (F 0 -N(d 0 ,v 0 )) 

(22) 


Equations (l 6) and (IT) are predictor equations (upper tilda) , (l8) and (19) are 
corrector equations, (20) to (22) are initial conditions, and N is a nonlinear 
algebraic operator^. The implementation procedure can be found in (ref. l). 

If frictional contact occurs, we need in addition a predictor equation for 
jl . Because of lack of space, this is not developed here. For the time being 
we adopt 



NUMERICAL RESULTS 


The analysis of an impact of a rectangular block on a rigid surface is 
performed (see .ref. 3, for comparison). Figure U shows the mesh. The data are 


density 

modulus of elasticity 
Poisson's coefficient 
dimensions 
time step 


p = 0.01 

E = 1,000 
V = 0.3 
L . W = 9 -9 
At = 0 . 002725 


f If variables are to be memorized at element integration points, as often in 
nonlinear problems, remember that these are moving when node 5 moves. 
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Newmark parameters 
(explicit predictor- 
corrector algorithm) 

initial velocity 

wave velocity 


Y = 0.5 0 = 0.25 


C d = {[E(1-V)] / [(1 + V)(1-2V). p]} 0 ' 5 = 366.9 


The time step is defined by the transit time for a dilatational wave to 
cross one element. The impact takes place at t = 0. Frictionless contact is as- 
sumed (M-^ = 0). This is introduced via isolated nodes, as shown on figure 4a. 

For the purpose of testing the new formulation, both node-to-node and distinct 
nodal positions are tested and yield the same results. 

The anticipated solution is shown on figure 4b. This exact solution has 
two constant zones separated by the dilatational wave front emanating from the 
initial impact. The circular wave front is a result of reflections off the free 
boundary . 

During the early steps of the computation, stresses in zone II are obtain- 
ed from the impulse equation applied to a one-dimensional situation ((J = c-p-v 0 ) 
Stress results shown on figure 5 -a confirm the validity of the new approach. 
Some overshoot appears, however, in the stress results of the lowest row of 
elements , probably due to the absence of a discrete impact condition in the 
algorithm. The deformed configuration at t = 0 . 0218 s . is shown on figure 5b. 


CONCLUDING REMARKS 


A new approach to contact problems involving friction in two-dimensional 
elastodynamics is proposed in this article. The basic idea obviously shows some 
analogies with local remeshing techniques, like the one proposed in (ref. 4). 

The treatment of friction via modified shape functions seems similar to the 
lines of thinking adopted in (ref. 5) for the treatment of shock waves. 

The proposed formulation is symmetric (both bodies are contactor and tar- 
get for each other) and satisfies compatibility of velocities over the contact 
interface (when possible), thus avoiding interpenetration. Completeness of the 
modified elements remains satisfied. Although no detailed comparison with dif- 
ferent approaches has been made as yet, the following advantages can be men- 
tioned : constant size of the system of equations (not true for local remeshing 
or Lagrange multiplier approach, important if an implicit solution is performed) 
and interface compatibility (not true in general for Lagrange multipliers or 
penalty methods ) . 

The extension of the method to several contacting nodes per element is 
possible, but it is not trivial. The extension to contact problems in three 
dimensional space and inclusion problems in 2-D is possible, but at the cost 
of losing interface compatibility. As already mentioned, nonlinear analysis may 
present minor difficulties because of the fact that integration points can move. 
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Further research is needed on the predictor algorithm for sliding with friction 
and impact-release conditions have to he added. 

Our main effort, at present, is directed towards testing the approach in 
problems involving friction. 
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Figure 1.- Contact problem. 
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Figure 2.- Modified shape functions 








